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This paper presents a canonical dual approach for solving a nonlinear population growth problem gov- 
I emed by the well-known logistic equation. Using the finite difference and least squares methods, the 

nonlinear differential equation is first formulated as a nonconvex optimization problem with unknown 
parameters. We then prove that by the canonical duality theory, this nonconvex problem is equivalent to a 
concave maximization problem over a convex feasible space, which can be solved easily to obtain global 
optimal solution to this challenging problem. Several illustrative examples are presented. 

■ 1- Problems and Motivations 

■ The logistic equation is a model of p opulation g r owth first published by Pierre Verhulst in 1838. The 
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continuous version of this model (see iMav et al.\ (119791) ) is described by the following first-order non- 
linear differential equation: 

f^.x(i-|,-c, a..i, 

where x(r) represents the numbers of individuals (population biomass) at time t, the real number r > 
M ' is the intrinsic growth rate of population increase, K is the carrying capacity, or the maximum number 

■ - -' of individuals that the environment can support, C > represents the constant harvesting rate. 

By using finite difference method, the discrete version of the logistic equation can be written as 

x,+i = x, + rx,(l-^)-C, + E„ (1..2) 

where is the process error, which is critically important to the associated dynamics equation; x^ 
is the stock abundance in year t. Due to the nonlinearity, this equation has remarkable non-trivial 
properties and therefore it must be handled with special methods. It is well-known that for certain given 
parameters, direct iterative methods for solving this nonlinear equation may lead to chaotic solutions. 
Such a problem was first studied on computer by the theoretical population biologist Robert May in the 
late 1960's. This equation and its variants still puzzle the mathematicians. 

In real-world appUcations, the stock biomass x cannot be observed directly, therefore, an observa- 
tional model can be introduced below 

I, =^x, =!, + £,, (1..3) 

(c) Institute of Mathematics and its Applications 2005; all rights reserved. 
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where I? is the estimation of abundance, I, is the abundance observed, q is the proportional coefficient, 
and £; is the observation error Base on this model, the least squares method for minimizing both the 
observation error e, and the process error Et leads to the following optimization problem. 



min EC-^X' + 7 E^' 

r=l ^ t=\ 

s.t. x,+ i =Xr + rxr(l - ^) - C, +Er, f = 1, . . . ,n - 1, 



(1..4) 
(1..5) 



where a is the penalty factor Let 
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the vector form optimization problem (^o) can be written as: 



a 



K 



x^Ax + (D - (1 + r)/?)x + C 



{^i) min Pi{ii,q,K,r) = \\qx^\\\ 

s.t. q > 0,K > 0,r> 0, 

xeM.",qeR,K eR,reR, 

where I e M", C e M""' are given vectors, A = {A^^.} {5^} G M«x(«-i)x'', and 

5'. = / ^ if'=;=^ 

I otherwise. 



By setting y ~x/K, can equivalent population dynamics model can be proposed in the following. 



K 



(^2) min P2{y,q,K,r)= qy - - +- ry^ Ay + {D - {I + r)R)y + 

K L 

s.t. 1 >y >0,^>0,/:>0,r>0, 

y GM",^eR,/:eR,r eM. 

Clearly, this fourth-order polynomial optimization problem with unknown parameters is in general non- 
convex, which could have many local minimizers. Due to the lacking of sufficient conditions for identi- 
fying global optimizers, traditional convex optimization theories and methods for solving this nonconvex 
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problem are very difficult. Actually, many nonconvex minimization problems in global optimization are 
considered to be NP-hard (see Gao and Sherali, 2009). 

Canonical duality theory is a potentially useful methodological concept which was developed orig- 
inally from complementary variational problems in nonconvex mechanics (see Gao and Strang, 1989). 
This theory is composed mainly of (1) a canonical dual transformation method, which can be used 
to reformulate the nonconvex primal problem as a perfect dual problem without duality gap; (2) a 
complementary-dual principle, which provides an analytical solution form to the primal problem; (3) 
a triality theory, which can be used to identify both global and local extrema. This theory has been 
used successfully for solving a large class of nonconvex/nonsmooth/discrete problems in computational 
biology, global optimization, phase transitions of solids, finite element methods for post-buckling of 
large deformed structures, and Euclidian distance geometry problems (see Gao et al, 2000 - 2012). 
The goal of this paper is to apply this theory for solving the nonconvex minimization problem {^2)- 
In the next section, we first introduce briefly the canonical duality theory with a simple example of a 
one dimensional double-well function optimization problem. Then we use the canonical dual transfor- 
mation to construct the canonical dual problem in Section 3. Furthermore, the form of the analytical 
solution is obtained from the criticality condition in Section 4. In Section 5, we present some numerical 
experiments. Concluding remarks are given in the last section. 



2. Canonical duality theory: A brief review 

The basic idea of the canonical duality theory can be demonstrated by solving the following general 
nonconvex problem (the primal problem (^) in short) 

(^): min jp(x) = -(x,Qx)-(x,f)+W(x)l, (2..6) 

xG,«fl |_ 2 J 

where Q e R"^" is a given symmetric indefinite matrix, f e M" is a given vector, (x,x*} denotes the 
bilinear form between x and its dual variable x*, W{x) is a general nonconvex function; and C M" 
is a given feasible space. 

The key step in the canonical dual transformation is to choose a nonlinear operator, 

^ = A (x) : X, -^SafZ W (2. .7) 

and a canonical function V : Sa^M. such that the nonconvex functional W(x) can be recast by adopting 
a canonical form W(x) = y(A(x)). Thus, the primal problem [iP) can be written in the following 
canonical form: 

(^): min {f(x) =y(A(x))-t/(x)}, (2..8) 

where ?7(x) = (x,f) — |(x,Qx). By the definition introduced in Gaol ( 2000l) . a differentiable function 
y(^) is said to be a canonical function on its domain if the duality mapping q = ^^(^) from <f„ to 
its range C W is invertible. Let (1^; 5) denote the bilinear form on K'^ Thus, for the given canonical 
function ), its Legendre conjugate V* [q) can be defined uniquely by the Legendre transformation 

y*(g) = sta{(^;g)-y(<^)| ^ G ^4, (2.-9) 

where the notation sta{g((^)| S, G So} stands for finding stationary point of g{^) on It is easy to 
prove that the following canonical duality relations hold on Sa x '^a- 
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g = VV(^) ^ ^=Vy*(g) ^V{^)+V*{q) = {^-q). (2..10) 

By this one-to-one canonical duality, the nonconvex term IV (x) = y(A(x)) in the problem {^) can be 
replaced by (A(x);5) — V*{q) such that the nonconvex function f (x) is reformulated as the so-called 
Gao and Strang total complementary function: 

S(x,g) - (A(x);g) -r(g)-t/(x). (2..11) 

By using this total complementary function, the canonical dual function P^{g) can be obtained as 

P^'ig) = sta{S(x,g) |xe JTJ 

= U^{g)-V*{q), (2..12) 

where (x) is defined by 

f/^(g)=sta{(A(x);g)-f/(x) I x€^a}- (2..13) 
In many applications, the geometrically nonlinear operator A (x) is usually quadratic function 

A(x) = i(x,D,x) + (x,b,), (2.. 14) 

where G M"^" and ^M."{k— 1, • • • ,p) are given. Let q = [qi,- ■ ■ ,gp]'^ . In this case, the canonical 
dual function can be written in the following form: 

P'iq) = -^{F{q),G-\q)F{q))-V*{q), (2..15) 

where G(g) = Q + I^^j and F{q) = f - I^i 5^^,. 
Let 

= {qe RPl G{q) y 0}. 
Therefore, the canonical dual problem is proposed as 

(^''): max{P''(g)| g £ ^+}, (2..16) 

which is a concave maximization problem over a convex set S^^^ C M.P. 

Theorem 2..1 (IgaoI (l2000h ) Problem {^'') is canonically dual to (^) in the sense that if g is a 
critical point of then 

x = G-Hg)F(g) (2..17) 

is a critical point of P{x) and 

P(x) = S(x,g)=P^(g). (2..18) 
If 5 € is a solution to {^'^), then x is a global minimizer of {^) and 

minP(x) = S(x,g)= max (2.. 19) 

"G^a qsy'+ 

Conversely, if x is a solution to (^), it must be in the form of ( I2..17l i for critical solution q of P^{q)- 
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To help explain the theory, we consider a simple nonconvex optimization in R": 

minP(x) = ^a{^\\xf-Xf-x^f, VxeM", (2..20) 

where a, A > are given parameters. The criticality condition VP(x) = leads to a nonlinear algebraic 
equation system in M" 

a(i||x||2-A)x = f. (2..21) 

Clearly, it is difficult to solve this nonlinear algebraic equation directly. Also traditional convex opti- 
mization theory can't be used to identify global minimizer However, by the canonical dual transforma- 
tion, this problem can be solved. To do so, we let ^=A{u) = i||xf -1 eR. Then, the nonconvex 
function W(x) = 5a(5 ||x|p — A)^ can be written in canonical form y(^) = jCtt^^. Its Legendre conju- 
gate is given by V* (g) = ja^^g^, which is strictly convex. Thus, the total complementary function for 
this nonconvex optimization problem is 

S(x,g) = (i||xf-A)g-ia-'g2_x^f. (2..22) 

For a fixed G M, the criticality condition VxS(x) — leads to 

gx-t=0. (2..23) 

For each g ^ Q, the equation ( 12. .23b gives x ~ f/g in vector form. Substituting this into the total 
complementary function E, the canonical dual function can be easily obtained as 

P'ig) = {S(x,g)|V,S(x,g) = 0} 

= -^-^a-'g^-Xg, Vg^O, (2..24) 

which has only one variable! The critical point of this canonical function is obtained by solving the 
following dual algebraic equation 

(a-ig + A)g2 = Vf. (2..25) 

For any given parameters a, X and the vector f G M", this cubic algebraic equation can be solved 
analytically to have at most three real roots satisfying gi ^ g2^ g?,, and each of these roots leads to 
a critical point of the nonconvex function P{x), i.e., x,- = f/g,-, ; = 1,2,3. By the fact that gi G ,5^^^ = 
{g eR \ g > 0}, then Theorem 1 tells us that xi is a global minimizer of P{x). Consider one dimension 
problem with a = I, X = 2, f = j, the primal function and canonical dual function are shown in Fig. 
[T] where, xi = 2.11491 is global minimizer of P{x), gi = 0.236417 is global maximizer of P^'ig), and 
P{xi) = -1.02951 (Seethe two black dots). 

The can onical duality theory was original developed for handling general nonconvex systems (see 
lOaol (I2OOOI) ). The canonical dual transformation can be used to convert a nonconvex problem into a 
canonical dual prob lem without duality gap, while the classical dual approaches may suffer from having 
a potential gap (see iRockafellan d 1 9991) ). The complementary-dual principle provides a unified form of 



analytical solutions ( I2..17l i to general nonconvex problems in continuous or discrete systems. 
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Fig. 1. Graphs of the primal function P(x) (blue) and its canonical dual function P^(q) (red). 



3. Canonical Dual Approach for Solving Nonconvex Problem (^2) 

In order to solve the nonconvex minimization problem (,^2) by the canonical duality theory, we need 
to fix the parameters r,K first and introduce a geometrical nonlinear measure Ai : R" ^M"-i defined 
by 

^ =Ai(y)-ry^Ay+(D-(l+r)/?)y + ^. 

The canonical function associated with this nonlinear measure is a quadratic function: Vi{^) — ^ajli^ll^ 
and the duality relation 

is invertible. Then the Legendre conjugate is defined by 

yr=max{^^g-y,(^)|^er'->} = ^g^g. 

On the same time, we rewrite the Unear inequality constraints < y, < 1 , / = 1 , • • • , n in the canonical 
form: 

yo(y-e) <0, 

where e = (1, . . . , 1) G R" is an one-vector, the notation sot := {siti,S2t2, ■ ■ ■ ,Sntn) denotes the Hadamard 
product for any two vectors s,t G R". Introduce a quadratic geometrical operator A2 : R" R" 

e = A2(x) =xo (x-e), 

and indicator function 



Viie) = 




ife <0, 
otherwise, 



then we know that, by theory of convex analysis, the subdifferential a G <3V2(e) is equivalent to the 
following KKT conditions: 

e <0, a^Q, e^a^Q. 
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Therefore, for the fixed parameter q,r,K, the inequahty constrained primal problem {^2) can be written 
in the following canonical form. 



I 









(^3) min |P3(y)=Vi(Ai(y))+y2(A2(y)) + 
Since the indicator function V2 is nonsmooth, its conjugate is defined by the Fenchel transformation 
y«(a) = sup{s-a - V2ie) \e e M«} = { « ^ J^^^ise. 
Thus, the primal function can be reformulated as the total complementary function: 



S(x,g,cT) = g^Mx)-Vng) + cj^A2{x)^vl{a) + \\qx^-f 

= V(g, a)G(g, (T)x(g, CT) - F^(g, (t)x - V,*ig) - V^{a) 



in which 



G(g,cj) = 2(rA^g + <?^// + Diag(a)), 
¥{q,a) = 2|l + (T-(D-(l+r)/?)^g, 

where // G R"^" is an identity matrix, and Diag (a) G R"^" denotes a diagonal matrix with {a,} (/ = 
1 , • • • ,n) as its diagonal entries. For a fixed g , (7, the criticality condition VxS (x, g , cr) =0 lead to the 
following canonical equilibrium equation: 

G(g,(j)x = F(g,(j), 

This equation can be solved analytically in the form of 

x = G-\g,a)F{g,a) (3..26) 

on the canonical dual feasible space defined by 

■ya = {{g,<j) e r("-"+"| cj ^ 0, detG(g,cj) /O}. 

Thus, the canonical dual problem can be formulated as 

(^3^) : sta|7^^'(g, a) - -^-F{g, afG-\g, a)F{g,a) - ^g\ + ^C^g + ^\^\ | g e ^„ 

Theorem 3.. 1 (Complementary-Dual Principle) For the given parameters q, K, and r, if (a, g) 
is a KKT point of {^^), then the vector 

y = G-i(g,a)F(g,a) (3..27) 
is aKKT point of (^3), and P3(y) = d). 
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Proof. Suppose that {c,,o) is a KKT point of {^f), by introducing Lagrange multiplier e to relax 
the inequality condition a ^ of the dual problem (.^^3 ), we have 

C 1 



ad I 





0, 


(3. .28) 




e, 


(3. .29) 


< 


0, 


(3. .30) 




0. 


(3..31) 


(1 


+ r)R. 





By the equation ( |3..29t and the KKT conditions (3..30I) and the complementarity condition j3..31l l. 
we have 

VP3(y)+2Diag(a)y-a = (3..32) 
0<y<l,a ^ (3. .33) 

^(y°y-y) = o, (3..34) 



where 



Therefore, 



y = (2(rA^g + <72// + Diag(a)))-i(2^| + a-B^g) 

is a KKT point of the primal problem {^3). 

Moreover, in terms of y = G^^{g, &)F{g, a), we have 

= -iF{m^G-'i-g,aMm-^f-g + ^C^-g + ^i^i 

= ix^(2(rA^g + ^2^ + Diag (d)))x -x^(2^1 + G-B^g) - ^fg + hz^g + 

= f{rfKy) + {Byf-g + U^g - ±f-g + f(q^H)y - 2{qyf^ + ^1^1 

= ^\\rf Ay + By +^\\'-+\\qy- ^\\'-=P,iy) 

This proves the theorem. □ 
For further discussion on extremality properties of the solution ( 13. .27b . we introduce the following 
feasible space 

= e ^„ I G(g,a) >- 0}. (3..35) 

Then, the extremality of the solution ( I3..27l i can be identified by the following theorem. 
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Theorem 3. .2 For the given parameters q,K,r, suppose that € is a critical point of the 

canonical dual function P^{g,&) and y = G^^{g, a)F{g, a). Then, y is a global minimizer of P3(y) on 
M" if and only if (g, a) is a global maximizer of P^{g,a) on S^^, i.e., 

P3(y)-minP3(y)^ max (g, a) = P3''(g, ff). (3..36) 

yeR" {g.a)€.y+ 

Remark 3..1 Whence the global optimal solution y of the nonconvex problem (^3) is obtained by 
the canonical dual approach, the optimal parameters q,r,K can be obtained by solving the following 
minimization problem 



min Pi[{q,K,r) 
s.t. q>Q,K>Q, r > 0. 



I 


' a 




+ 2 



2 



By combining together problems (^3) and {^4), the global optimal solution {y,q,r,K) of the problem 
(^2) can be obtained by certain alternative iteration method. Then, the global minimizer to the problem 
(^1) isx=^y and we havefi(x,^,^,r) ^ K^P2{y,q,K,r). 

4. Applications 

We now list a few examples to illustrate the applications of the theory presented in this paper. 



Ex ample 1. The data we used in Table 1 is from Quinn T and Deriso R.B.'s book (see lOuinn & Deriso 
(ll999h V 



By the canonical duality method presented in this paper, the optimal parameters we obtained are 
q = 0.01 12, K 14777.4, r = 0.1875, and the optimal solution y to model (^2) is 

y = {0.87804,0.595806,0.465162,0.407183,0.368356,0.353525,0.353367,0.366521,0.368154, 
0.372081,0.360713,0.391465,0.409868,0.432046,0.445403,0.458043,0.467348,0.477866, 
0.490154,0.509345,0.542681,0.476886,0.43092,0.386888,0.351872,0.296357,0.245213}. 

The optimal objective function value is therefore P2{y,q,K,r) = 0.0698. Figure 2 presents a clear view 
of the tendency during these 27 years. 

Example 2 The data w e used in this example is the New Zealand rock lobster during 46 years given 



m 



Polacheck ef g/.l (119931) 



By the canonical duality method, the optimal parameters for the given data are <7 = 0.0005, K = 
28238.1, r = 0.5729, and the optimal solution y to model {^2) is 

y = {0.539323,0.514615,0.515001,0.537057,0.548343,0.566006,0.542414,0.557294, 
0.586471,0.591127,0.535325,0.496853,0.40859,0.403013,0.413048,0.429969, 
0.473993 , 0.49995 , 0.488957, 0.490067, 0.498439, 0.484909, 0.450408 , 0.464355 , 
0.448175,0.447358,0.432635,0.401602,0.438707,0.429822,0.411511,0.461489, 
0.468958,0.494921,0.517036,0.497474,0.480722,0.491527,0.489311,0.49443, 
0.459034, 0.40785 , 0.368474, 0.275356, 0.275089, 0. 166475} . 

The optimal objective function value P2 (y, q, K, r) ~ 2. 1596 (see Figure|3]l. 
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Fig. 2. Stock abundance during 27 years. 



Table 1 Yield (Catch) and catch-per- 
unit-effort I ( CPUE) for each of year 



Year C (Catch) I (CPUE) 



1 


3706 


82 


2 


2662 


59 


3 


2055 


46 


4 


1658 


37 


5 


1379 


31 


6 


1171 


26 


7 


1011 


22 


8 


884 


20 


9 


781 


17 


10 


696 


15 


11 


85 


17 


12 


111 


22 


13 


142 


28 


14 


178 


36 


15 


215 


43 


16 


253 


51 


17 


289 


58 


18 


321 


64 


19 


349 


70 


20 


371 


74 


21 


1451 


73 


22 


1341 


67 


23 


1263 


63 


24 


1206 


60 


25 


1162 


58 


26 


1129 


56 


27 


1103 


55 
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Fig. 3. Stock abundance during 46 years. 



5. Concluding remarks 

We have presented a canonical dual approach for solving a challenging population growth problem in 
discrete dynamical systems. By using the finite difference and the least squares methods, this nonlinear 
problem is reformulated as a nonconvex optimization problem with inequality constraints. The problem 
has extensive applications in ecology, neural networks, medicine, economics, statistical physics, and 



Table 2 Catch and catch-rate data 



Year 


Catch 


CPUE 


Year 


Catch 


CPUE 


1945 


809 


3.49 


1968 


4975 


1.53 


1946 


854 


3.38 


1969 


4786 


1.32 


1947 


919 


3.18 


1970 


4699 


1.45 


1948 


1360 


3.56 


1971 


4478 


1.40 


1949 


1872 


1.79 


1972 


3495 


1.09 


1950 


2672 


4.35 


1973 


3784 


1.23 


1951 


2834 


2.33 


1974 


3643 


1.12 


1952 


3324 


2.57 


1975 


2987 


0.92 


1953 


4160 


2.88 


1976 


3311 


1.02 


1954 


5541 


3.85 


1977 


3237 


1.0 


1955 


5909 


4.16 


1978 


3418 


1.05 


1956 


6547 


4.34 


1979 


4050 


1.09 


1957 


5049 


3.70 


1980 


4190 


1.02 


1958 


4447 


2.37 


1981 


4058 


1.01 


1959 


4018 


2.46 


1982 


4331 


0.98 


1960 


3762 


2.06 


1983 


4385 


1.01 


1961 


4042 


2.21 


1984 


4911 


0.85 


1962 


4583 


2.19 


1985 


4856 


0.84 


1963 


4554 


2.44 


1986 


4657 


0.81 


1964 


4597 


2.14 


1987 


4500 


0.84 


1965 


4984 


2.18 


1988 


3128 


0.68 


1966 


5295 


2.13 


1989 


3318 


0.62 


1967 


4782 


1.86 


1990 


2770 


0.54 
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much more. Due to the nonconvexity of the target function with unknown parameters, the problem 
is fundamentally difficult and can't be solved by traditional direct methods. However, by the canon- 
ical duality theory, global optimal solution is obtained for the first time to this well-known problem 
over the whole time domain. Applications illustrated that this theory is efficient for solving large-scale 
population growth problems. 
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